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The equation of state (EOS) for nuclear and neutron rich matter is investigated in a Relativistic 
Mean Field (RMF) model. New shell states are found that minimize the free energy per baryon, 
calculated in a spherical Wigner-Seitz (WS) approximation, over a significant range of baryon den- 
sities. These shell states, that have both inside and outside surfaces, minimize the Coulomb energy 
of large proton number configurations at the expense of a larger surface energy. This is related to 
a possible depression in the central density of super heavy nuclei. As the baryon density increases, 
we find the system changes from normal nuclei, to shell states, and then to uniform matter. 
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I. INTRODUCTION 

The equation of state (EOS) for dense nuclear matter 
plays an important role in understanding the structure of 
neutron stars and for simulations of core collapse super- 
novae. Neutron stars probe the zero temperature EOS 
while supernovae involve the EOS at temperatures up to 
tens of MeV. Below the neutron drip density, « 2 x 10^^ 
g/cm.^^ where free neutrons appear, knowledge of the EOS 
can be extrapolated from the properties of laboratory 
atomic nuclei. In contrast, current terrestrial experi- 
ments only constrain the EOS above neutron drip density 
with rather big uncertainties. Recently, there have been 
efforts from heavy ion collision experiments [l|, X-ray 
observations of isolated neutron stars [3] , and theoretical 
calculations to understand the EOS beyond the neutron- 
drip point. 

There have been many calculations of nonuniform mat- 
ter at sub-nuclear density, below 3 x 10^"' g/cm^. Negele 
and Vautherin, using a Skyrme Hartree-Fock calculation 
with two-body potentials Q, obtained the ground state 
within the spherical Wigner-Seitz (WS) approximation. 
In this approximation, the unit cell of a crystal lattice is 
modeled as a sphere. They found, as the system became 
more neutron rich and the density increased, neutrons 
escaped and the system approached a uniform state near 
nuclear density. Since then there have been many in- 
vestigations Q with more sophisticated interactions and 
more complicated lattice configurations, such as cylin- 
ders or plates. These non-spherical 'pasta' phases, which 
seem to appear within a significant range of sub-nuclear 
densities are relevant for the structure of neutron star 
crusts d, @| and the dynamics of supernovae [T] . 

For hot nuclear matter, assuming that the spherical 
WS cell still remains a good approximation to the real- 
istic lattice structure, people have studied the EOS with 
different models. Using a phenomenological compress- 
ible liquid-drop model, Lattimer and Swesty 0] produced 
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an equation of state for hot dense matter that has been 
widely used in supernova simulations. Later, H. Shen 
et. al. constructed an equation of state based on 
Thomas-Fermi and variational approximations to a rel- 
ativistic mean field energy functional. Neither method 
takes into account the shell structure of finite nuclei or 
explores the full range of density distributions possible 
even in the spherical WS approximation. 

In recent years, relativistic mean field models (RMF) 
have provided a consistent description for the ground 
state properties of finite nuclei, both alon g an d far away 
from the valley of beta stability flO, [O, PJ . These 

models incorporate the spin-orbit splitting naturally and 
the relativistic formalism provides a framework to extrap- 
olate the properties of non-uniform and uniform nuclear 
matter to high densities. Furthermore there is a close re- 
lation between phenomenalogical RMF models, that sim- 
ply fit parameters to properties of finite nuclei, and more 
systematic effective field theory approaches that enumer- 
ate all of the possible interactions that are allowed by 
symmetries. 

In this work, we develop a RMF code to explore the 
equation of state for non-uniform and uniform nuclear 
matter at finite temperatures with a range of proton 
fractions and wide range of baryon densities. For non- 
uniform nuclear matter, we use the spherical WS approx- 
imation to find the stable state for a charge neutral sys- 
tem of neutrons, protons and electrons. In this paper, 
we will focus on the matter at a low temperature T =1 
MeV. It is easier for us to find the ground state at T = 
by extrapolating the theory at a low temperature T ^ 
to T = 0, because the occupations of nucleon levels vary 
more smoothly at finite temperature. In the future we 
will also study the EOS at higher temperatures. 

For nonuniform matter we find new shell states which 
minimize the free energy per baryon over a significant 
density range. Shell states have inside and outside sur- 
faces and they can minimize the Coulomb energy of high 
Z (large proton number) configurations at the expense of 
a larger surface energy. These shell states may be related 
to the appearance of a central depression in the density of 
super heavy nuclei [l^ because of their large coulomb en- 
ergies. The appearance of shell states may significantly 
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change transport properties such as the shear viscosity 
and shear modulus of neutron rich matter. 

We compare the transition density, that we find in 
nonuniform RMF calculations, to that predicted from a 
stability analysis of uniform nuclear matter. We start 
from uniform matter in equilibrium, and then find the 
critical density when the longitudinal collective mode, 
that describes density oscillations, becomes unstable. 

The paper is organized as follows. First in Section [Til 
we give a brief description of the RMF theory at finite 
temperature and explain the WS approximation for non- 
uniform matter. Next, Section lll BI presents the collective 
mode analysis for uniform matter in both the Vlasov for- 
malism and using the RPA approximation. In Section 
mil we present results for the EOS of nonuniform mat- 
ter. Finally in Section [TVl we summarize our results and 
comment on future work. 



II. FORMALISM 

The RMF theory and its applications have been re- 
viewed in a number of places, see for example Refs. [12, 
[Tsl . nil. The basic ansatz of the RMF theory is a La- 
grangian density where nucleons interact via the ex- 
change of sigma- (cr), omega- (lo), and rho- (p) mesons, 
and also photons (A). 
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Here the field tensors of the vector mesons and the elec- 
tromagnetic field take the following forms: 
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In charge neutral nuclear matter composed of neu- 
trons, n, protons, p, and electrons, e, there are equal 
numbers of electrons and protons. Electrons can be 
treated as a uniform Fermi gas at high densities. They 
contribute to the Coulomb energy of the npe matter and 
serve as one source of the Coulomb potential. 

The variational principle leads to the following equa- 
tions of motion 

[a-p + y(r) + /3(m + S'(r))]V'i = eiV^i (3) 

for the nucleon spinors, where 



S{y) = g^a 



(4) 
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for the mesons and photons, where the electrons are in- 
cluded as a source of Coulomb potential. The nucleon 
spinors provide the relevant source terms: 

Ps = J2^ 

At zero temperature, the summations run over the va- 
lence nucleons only, since rii = Q^ep — £i) where ei? is 
the Fermi energy of the nucleons. At finite temperature, 
Fermi-Dirac statistics prescribes the occupations of pro- 
tons and neutrons as follows: 
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where p, is the chemical potential for neutron (proton). 
In practice, we include all levels with Ui > 10~^. 

Since the systems under consideration have tempera- 
tures of, at most, tens of MeV, we neglect the contribu- 
tion of negative energy states, i.e., the so-called no sea 
approximation. In a spherical nucleus, there are no cur- 
rents in the nucleus and the spatial vector components 
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and A,, vanish. One is left with the time- 



like components. Charge conservation guarantees that 
only the 3-component of the isovector po,3 survives. The 
above non-linear equations are solved by iteration within 
the context of the mean field approximation whereby the 
meson field operators are replaced by their expectation 
values. 



A. Nonuniform Matter with Wigner-Seitz 
approximation 

The spherical Wigner-Seitz (WS) approximation is 
used to describe non-uniform matter. In general the unit 
cell of a crystal lattice is a complex close-packed polyhe- 
dron. This is approximated by a spherical cell of the same 
volume. We apply boundary conditions on the wave func- 
tions at the edge of WS unit cell. To achieve a uniform 
density distribution for a free neutron gas, we require 
that at the cell radius, all wave functions of even par- 
ity vanish, and the radial derivative of odd parity wave 
functions also vanishes Q . 

In the WS approximation, the lattice Coulomb en- 
ergy consists of contributions from neighboring unit cells. 
This correction to the WS Coulomb energy is important 
for determining the stable configuration of WS cells and 
the transition density to uniform matter, especially when 
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the system has a large proton fraction. We include the ex- 
act Coulomb energy in calculating the free energy of WS 
cell of radius -Rc- Following the treatment in Ref. p^fl^. 
we calculate the Coulomb energy per unit cell as, 
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We also include the contribution of free electrons. Ex- 
plicit formulas for the electron contribution to the free 
energy per baryon /e, and entropy Se of the electron gas 
are presented in Ref. [§]. Finally, the complete free en- 
ergy per baryon is 



/ = + fe- 
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where a is lattice constant defined by — Vceii — ^R^, 
and 



G = hA + kB + IC. 



(8) 



Here h, k, and I are integers, and A, B, and C are the 
primitive transformation vectors of the reciprocal lattice. 
The prime on the sum means that the point G = is 
excluded. The form factor I^ki is given by 
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Oyamatsu et al. [l6| assumed a uniform proton density 
inside the nucleus and found that the stable configuration 
is a Body-Centered Cubic (BCC) lattice. Using realistic 
proton density distribution, we also find that a BCC lat- 
tice gives the lowest Coulomb energy. In principle, one 
needs to worry about the screening of Coulomb interac- 
tions by electrons. It turns out that screening contribu- 
tions to the Coulomb energy are very small [5| and we 
will neglect them in our calculations. 

With the approximations specified above, it is conve- 
nient to perform a self-consistent relativistic mean field 
Hartree calculation for nuclear wave functions inside a 
WS cell of radius Rc, for given average baryon density 
Pb, proton fraction Yp and temperature T. The nucleon 
number inside the WS cell is therefore A = AttR^pb/S 
and proton number Z = YpA. The total energy of the 
WS cell, including the exact lattice Coulomb energy, is 
obtained as follows, 

Eb = Enucleon + E„ + Ep + + Wc - TuA 

= Y^Cin., - e j ppAo{r)(fr ~ i j (fr{g„aps{r) 
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(10) 

The nucleon contribution to the entropy is given by the 
usual formula, 

Sb = -kB^[n.Mn.,) -f (1 - 7i,)ln(l ~ nO] , (11) 

i 

where is given in Eq. ([6]). With Eqs. (fTO]) and pT|) . 

it is easy to obtain the nucleon contribution to the free 
energy per baryon fb, 



h = Fb/A = (Eb - TSb)/A. 



(12) 



B. Stability analysis 

1. Semiclassical collective modes analysis 

The unstable collective modes of uniform npe matter 
at finite temperature can be investigated in the RMF 
by using a Vlasov formalism [l3|. Here we give a short 
summary to make this paper more self-contained. 

The distributions of particles (antiparticles) denoted 
as -|- (-), at position r, time t, and momentum k (— k) 
are described via the phase space distribution functions 
as following: 



/(r,k,i). 



diag(/p±,/„±,/e±)- 
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The time evolution of distribution functions are deter- 
mined by the Vlasov equation: 



df, 



at 



+ {fi±,K±} = 0, i =p, n, e. 



(15) 



where { , } denotes the Poisson brackets. hi± is the one 
body Hamiltonian, which can be derived from Eq. ([T]) 
for neutrons and protons. For electrons, h is the nor- 
mal QED Hamiltonian. Small deviations Sf± of the dis- 
tribution functions around the equilibrium state can be 
calculated with generating functions 



S{r,k,t)± = diag(S'p±,S'„±,S'e±), 



such that 



sf± = {^±,/o±} - {s±,e} 



dfo± 
dk^ ' 



(16) 



(17) 



where /o± are equilibrium distribution functions. In 
terms of the generating functions, one can obtain lin- 
earized Vlasov equation for 8f±. From the equations of 
motion for meson and photon fields, one can also get the 
linearized equations for the oscillation fields. 

The dominant collective modes at low temperatures 
are the longitudinal modes (with momentum p and fre- 
quency 17), which can be described by the ansatz 




^^±(fc,cos0) 



(18) 



where j — p, n, e, — oj^, p^, A^ represents the vector 
fields and photon fields. 6 is the angle between p and k. 
After transforming the unknown variables to the density 
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oscillations of p, n, e, we can obtain the following matrix 
equation for the eigenmodes 



D{n) 



Spnp 



0, 



(19) 



where Spf^^ and Sp^^ are the amplitudes of the oscillating 
scalar densities of protons and neutrons, respectively, and 
^Pnp, ^PVLni and ^poe are the amplitudes of the oscillating 
proton, neutron, and electron densities. The entries of 
the matrix D(f2) are given in Ref. 
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The spectrum of collective modes, or the dispersion 
relation f7 = ri(p), is determined by the eigen condition 
D(yt) — 0. The instability of collective modes can be 
deduced directly from the static limit of the dispersion 
relation, i.e., if D{Q. = 0) < 0, the collective modes will 
be unstable and exponentially growing. This indicates a 
phase transition to non-uniform matter. 



2. RPA at finite temperature 

We now perform a second very similar stability analy- 
sis using the random phase approximation (RPA). This 
provides an independent check of the Vlasov approach. A 
detailed account of the finite temperature RPA method 
for Quantum Hadrodynamics can be found in Ref. p^ . 
Again we present some information to demonstrate the 
idea. The finite temperature Feynman rules yield the 
(s, s') thermal component of the causal polarization in- 
sertion for particle species i {i — n,p, e) as 
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^(-l)(^+^')Tr[G(^-;)(fc)J,(^) 

'G^i;:hk+q)j.{i)], 



(20) 



where G^^'^^-^ is the (s, s') thermal component of Hartree 

propagators of species i. For 7(photon), w, and p meson, 
= and the polarization is a tensor n^^^. For the 

a meson = 1, and the polarization is a scalar H^,. If 
— 7p and — \ (or vice versa), the polarization is 

mixed 11^^^^. 

At finite temperature, Dyson's equation for the proper 
longitudinal causal polarization is a matrix with particle, 
Lorentz, and thermal indices. It can be written as 



where A p (k) is the normal Feynman propagator and 
is the thermal matrix defined as 



cosh 9 — sinh 9 
sinh 9 cosh 9 



(23) 



The angle 9 is given by sinhfl = [exp/3|fco| — 1]~^/^, with 
inverse temperature /?. The thermal vector meson prop- 
agator has the same structure. 

At finite temperature we can also define 



(24) 

which allows us to derive the dielectric function by the 
real part of its (11) component. 
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detRe(£i((7^)") = det{l - A^Ren^^"^} 
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where we have used Rell^*-^^'' = Rell^. The poles of the 
dielectric function define the collective mode excitations 
of uniform matter. The dielectric function in the static 
limit determines the stability of the modes. When 



£L(go -0,g) <0, 



(26) 



the system is unstable against small amplitude longitu- 
dinal collective oscillations (density fluctuations). This 
indicates a phase transition to non-uniform matter. 

The above equation is for one particle species and one 
boson interaction. This can be generalized to several 
particle species with multiple interactions. For a system 
consisting of neutrons, protons, and electrons, the inter- 
action matrix and polarization matrix are given as 
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The explicit formula for these polarizations are listed in 
Ref. [is] . For the Lagrangian in Eq. ([1]) , the correspond- 
ing meson propagators are 

47ra 
.,2 



9a 



Si = m 



(21) 



where is the lowest order thermal matrix for the me- 
son/photon propagator. For example, for scalar meson 
propagator is 



q^ 



T*2 ■ 



where the effective masses are defined as 



(29) 
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I- 2g2(J + 3g3<J , 
icstol (30) 
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III. RESULTS 

In this work we use the NL3 effective interaction [2^ 
that has been successful in reproducing ground state 
properties of stable nuclei and the saturation properties 
of symmetric nuclear matter. The values of parameters 
for the NL3 effective interaction are listed in Table HI In 
this section, we present results in the spherical WS ap- 
proximation for non-uniform matter at a temperature of 
T = 1 MeV. 

To determine the minimum free energy per baryon / 
in Eq. (|13p at a specified baryon density ps, one must 
search over the cell radius Re- When Rc is large and pB 
is beyond the neutron drip density, one will need to take 
into account a large number of levels, since the nucleon 
number is ^ = AttR'^pb /•i- The fact that the matter is 
at finite temperature will require one to consider even 
more levels. For example, when pB — 0.080 fm~'^, Yp 
= 0.4, and Rc = 23.5 fm, there are 4892 nucleons inside 
one WS cell and one needs to include 427 (419) neutron 
(proton) levels. It is hard to achieve self-consistency for 
the mean fields with a large number of levels. To ensure 
the convergence of the self-consistent iterations for the 
mean fields, it is important to have a good initial guess for 
the mean field potentials. Here our strategy is to employ 
the convergent potentials for a nucleus with a smaller Rc 
as the starting guess for a nucleus with a slightly larger 
Rc- In this way, new shell states with large cell radii 
are found. These can minimize the free energy over a 
significant density range. 



per baryon become degenerate, as shown in upper right 
panel of Fig. [TJ Below 0.049 fm~^, the first minimum 
with smaller cell radius corresponding to a normal nu- 
cleus is the absolute minimum and therefore the equilib- 
rium state {e.g., panel a of Fig. [T]). On the other hand, 
above 0.049 fm~"^, the second minimum corresponding 
to a shell state is the true equilibrium state {e.g., panels 
c and d of Fig. [T]). Therefore 0.049 fm"'^ is the critical 
baryon density when the density distribution inside the 
WS cell changes from a normal nucleus to a shell state. 
In the lower right panel of Fig. [21 we also show the uni- 
form neutron and proton density distributions by dotted 
lines when pB = 0.090 fm~'^. As will be shown later, 
0.090 fm^'^ is the transition density from a shell state to 
uniform matter. 
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TABLE I: NL3 effective interaction. The nucleon masses are 
M — 939 MeV for botli protons and neutrons and C3 = 0. 

Qa Op 92 gs rria rUp 

(fm~^) (MeV) (MeV) (MeV) 

10.217 12.868 4.474 -10.431 -28.885 508.194 782.5 763 

We show in Fig. [T] an example of how shell states were 
found. Here the free energy per baryon versus WS cell 
radius is displayed for baryon densities of 0.040, 0.049, 
0.080, and 0.090 fm~^ when Yp =0.4. For each baryon 
density, there are two minima of free energy per baryon 
at different WS radii, denoted as open circles. The neu- 
tron and proton density distributions of the WS cells with 
these two radii are shown in Fig. [2l The first minimum 
at smaller WS cell radius corresponds to a nucleus with 
a normal density distribution. The second minimum at 
larger WS cell radius corresponds to a shell shaped den- 
sity distribution. Here the nucleus has both outside and 
inside surfaces and the neutron and proton densities are 
only non-vanishing at intermediate r. We call the first 
minimum a nucleus state and the second minimum a shell 
state. For the shell state, the sum of neutron and proton 
densities (at intermediate r) is close to the saturation 
density of nuclear matter ~ 0.15 fm~'^, as one expects 
from nuclear saturation. 

At pb = 0.049 fm^'', the two minima of free energy 



FIG. 1: (Color online) Free energy per baryon versus Wigner- 
Seitz cell radius Rc for baryon densities of ps — 0.040 (a), 
0.049 (b), 0.080 (c), and 0.090 fm"^ (d). The proton fraction 
is Yp = 0.4 and the temperature is T = 1 MeV. The two 
circles denote minima for the nucleus (small Rc) and shell 
states (larger Rc). 

There is an intuitive interpretation of shell states by 
considering the competition between the surface and 
Coulomb energies. When the Coulomb energy is more 
than half of the surface energy, the nucleus can decrease 
its Coulomb energy, via increasing its average radius and 
surface area, so that it minimizes the total energy. One 
possibility is to push the nucleons out spherically and 
form a shell shape. Note, this change is possible even 
if one assumes spherical symmetry. The nucleus could 
also deform its shape to a non-spherical configuration to 
minimize the total energy. It may be that, for very neu- 
tron rich matter in certain regions of density, the most 
stable configuration is non-spherical [2l[ . Therefore shell 
states may be related to the appearance of complex pasta 
phases [H, [21I, . Perhaps one can think of shell states 
as "spherical pasta" . 

To get a more quantitative idea for the appearance of 
shell states, we consider a simple model. Suppose there 
is a spherical WS cell with average baryon density pB = 
0.071 fm^'^ and Yr, = 0.5. Assume the nucleons have a 
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FIG. 2: (Color online) The neutron (solid line) and proton 
(dash line) density distributions for the Wigner-Seitz cells 
with baryon densities of ps = 0.040 (a), 0.049 (b), 0.080 (c), 
and 0.090 fm"^ (d). The proton fraction is Yp = 0.4 and T 
— 1 MeV. The red curves give density distributions for shell 
states. The black curves give density distributions for nucleus 
states. The corresponding cell radii are shown in each panel. 
The dotted lines in the lower right panel are neutron (upper) 
and proton (lower) densities for uniform matter. 



uniform density equal to the saturation nuclear density. 
There can be either a uniform density nucleus or a uni- 
form density shell state with a hole in the center. The 
schematic model is shown in the upper panel of Fig. [3l 
The relevant contributions to the energy are the surface 
and Coulomb energies. For each fixed WS cell radius, the 
radius of the outside, and perhaps inside surface of the 
nucleus can be adjusted to minimize the total energy. By 
this procedure of minimization, we can get the Coulomb 
energy, surface energy, and total energy of the cell for a 
range of cell radii as shown in the lower panel of Fig. [3l 
When the cell radius is small, the configuration minimiz- 
ing the total energy is that of a normal nucleus. When 
the cell radius increases, so that the Coulomb energy is 
larger than half of the surface energy, the nucleons will 
be pushed out from the center by Coulomb repulsion and 
form a shell. This transition decreases the Coulomb en- 
ergy while the surface energy is increased. For large WS 
radii, the change in Coulomb energy dominates over the 
change in surface energy and the total energy of the cell 
will decrease. Therefore the appearance of a shell state 
reduces the Coulomb energy and minimizes the total en- 
ergy at the expense of a larger surface energy. Now there 
are two minima in the total energy as a function of WS 
cell radius. The nucleus state has a radius of 12.52 fm 
and has "^Z =:^^'*292. The shell state has a radius of 
20.08 fm and ^Z^ 2408^^204. Here the shell state is more 
stable since it has a lower energy. 

There may be a similar situation for super heavy nu- 
clei. For the predicted doubly magic nucleus ^^^120, dif- 
ferent mean field model calculations show central depres- 
sions in the nuclear density fl^. The shell state we find 
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FIG. 3: (Color online) Upper panel: nucleon density distri- 
butions for nucleus (solid black curve) and shell state (solid 
red curve) in a schematic model for WS cell with ps = 0.071 
fm~'^ and Yp = 0.5; dotted line is electron density. Lower 
panel: the total energy per baryon Et/A, surface energy per 
baryon Es/A and Coulomb energy per baryon E^/A versus 
WS cell radius, with the same conditions as above. 



can be considered as an even larger super heavy nucleus 
with a central depression: the Coulomb repulsion is so 
large that the protons are pushed out from the center 
completely and locate in the outer region of the cell. Note 
that the neutrons then follow the protons in order to re- 
duce the symmetry energy. 

In Fig. m the atomic A and proton Z numbers for nu- 
cleus and shell states (upper panel) , and the correspond- 
ing WS cell radii (lower panel) are shown as functions of 
the baryon density. The proton fraction is Yp = 0.4. At 
low baryon densities, the WS cell radius of the shell state 
is about twice that of the nucleus state, e.g., 26 fm vs. 
14 fm at pb = 0.040 fm"-^. The cell radii of shell states 
decrease with increasing density, while those of nucleus 
state increase with density. They approach one another 
at high baryon densities: 24 fm vs. 22 fm at 0.090 fm^'^. 
For both nucleus and shell states, A and Z increase with 
increasing baryon density. 

In the upper panel of Fig.[5l the free energy per baryon 
for non-uniform matter (both nucleus state and shell 
state) and uniform matter are shown as function of den- 
sities with Yp = 0.4. In the lower panel, the free energy 
difference between the nucleus and shell states is shown 
for various baryon densities. This difference is of order 
tens of KeV and is maximum around 0.070 fm~'^. When 
the baryon density pB < 0.049 fm^^^, the nucleus state 
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is the most stable. Between densities of 0.049 to 0.090 
fm~'^, the shell state becomes slightly more stable than 
the nucleus state. When ps > 0.090 fm^'^, the matter 
becomes uniform. 




FIG. 4: (Color online) The atomic and proton numbers (panel 
a), and the WS cell radii for nucleus (squares) and shell states 
(filled circles) (panel b) shown as functions of the baryon den- 
sity for Yp = 0.4. 
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FIG. 5: (Color online) Upper panel (a): the free energy per 
baryon for non-uniform matter, nucleus (solid line) and shell 
states (dot-dash line), and uniform matter (dashed line) are 
shown as function of the baryon density for = 0.4 and T 
= 1 MeV. Lower panel (b): free energy difference per baryon 
between nucleus and shell states versus baryon density. 

The above analysis for Yp = 0.4 can be generalized for 
different proton fractions Yp. For Yp = 0.5, the minimiza- 
tion of free energy per baryon over the WS cell radii with 
baryon density of 0.045, 0.053, 0.073, and 0.090 fm^^ ^^^.^ 
shown in Fig.[Sl Again there are two minima correspond- 
ing to nucleus state and shell state respectively. The 
neutron and proton density distributions associated with 
the two states are shown in Fig.[71 The transition density 
from nucleus state to shell state is 0.053 fm^^. In Fig. [51 



the atomic and proton numbers (upper panel), and the 
WS cell radii (lower panel) for nucleus state and shell 
state are shown for various baryon densities. In Fig. [51 
the free energy per baryon for non-uniform matter (both 
nucleus state and shell state) and uniform matter are 
shown as function of densities. When the baryon den- 
sity is smaller than 0.053 fm~^, the nucleus state is the 
most favorable state. Between 0.053 to 0.090 fm~'^, the 
shell state become slightly more stable than the nucleus 
state. When the density is above 0.090 fm~^ (same as 
that when Yp — 0.4), the uniform matter is the stable 
state. 

Similar calculations were carried out for Yp = 0.3 and 
the results are presented in Figs. [TOl HH [121 and [T31 We 
note that the transition density to the shell state is 0.043 
fm~'^ and the transition density to uniform matter is 
0.084 fm~^. For Yp = 0.2, the results are presented in 
Figs, m [TH [TH and [m The transition density to shell 
state is 0.044 fm~^ and to uniform matter is 0.071 fm~^. 
The matter is very neutron rich at Yp — 0.2 and the nu- 
cleon density distributions of shell states have a distinct 
feature, compared to results for Yp = 0.3, 0.4, and 0.5. 
The neutron densities are not vanishing in the center and 
on the edge of the shell nuclei, as shown in Fig. [151 The 
neutron chemical potential is greater than zero. There- 
fore there are appreciable numbers of free neutrons inside 
the cell. 

Comparison between Fig. [9l and Fig. [TTl indicates that 
shell states will start to appear at smaller density for 
smaller Yp. Shell states involve a competition between 
surface and Coulomb energies. Neutron rich matter, with 
only a few protons, has a small binding energy and there- 
fore a small surface energy. This small surface energy 
makes it easier to form shell states. Furthermore, we no- 
tice that at smaller Yp the transition to uniform matter 
appears at smaller density. So in fact the density range 
for the shell state has decreased for smaller Yp. 

Up to now, we have examined the EOS in the WS ap- 
proximation for nonuniform and uniform nuclear matter 
with Yp = 0.2, 0.3, 0.4, and 0.5. A universal picture 
of phase transitions from normal nuclei to shell nuclei, 
then to uniform matter with increasing baryon densities 
emerges. It is beneficial to quantitatively compare the 
results to other analytical methods. As presented in Sec- 
tion. IIIBl using the eigencondition £)($! = 0) = from 
the semiclassical collective mode analysis and Eq. 
from the finite temperature RPA analysis respectively, 
one can obtain the phase transition density to uniform 
matter. In Fig. 1181 the phase transition densities to uni- 
form matter at T = 1 MeV are shown for different pro- 
ton fractions. The transition densities obtained from the 
above two analyses are plotted in comparison with the 
Hartree WS calculations. At Yp = 0.3, these two anal- 
yses give the transition density close to the WS approx- 
imation. When Yp > 0.3, the two analyses give slightly 
smaller transition densities. When Yp < 0.3, the two 
analyses give bigger transition densities. The latter dif- 
ference may suggest existence of non-spherical configura- 
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tions in the region between the spherical WS transition 
density and the collective modes/RPA transition density 
for non-uniform neutron rich matter. 
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FIG. 6: (Color online) Free energy per baryon versus Wigner- 
Seitz cell radius for baryon densities of pB = 0.045, 0.053, 
0.073, and 0.090 fm"^, proton fraction Yp = 0.5 and T = 1 
MeV. The two open circles denote minima corresponding to 
the nucleus (left) and shell (right) states. 
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FIG. 7: (Color online) The neutron (solid line) and proton 
(dash line) density distributions for the Wigner-Seitz cells 
with baryon density ps = 0.045, 0.053, 0.073, and 0.090 fm"^, 
with proton fraction Yp — 0.5 and T = 1 MeV. The red curves 
give density dist. of shell state. The black curves give density 
dist. of nucleus state. The corresponding cell radii are shown 
in each panel. The dotted line in the lower right panel is the 
neutron or proton density of uniform matter. 



IV. SUMMARY 

In this paper 'we applied a Relativistic Mean Field 
model to study the equation of state for non-uniform and 
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FIG. 8: (Color online) The upper panel (a) shows atomic 
(solid line) and proton (dashed line) numbers, while the lower 
panel (b) shows WS cell radii, for nucleus (squares) and shell 
states (filled circles) as functions of the baryon density for Yp 
= 0.5. 
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FIG. 9: (Color online) Upper panel (a): the free energy per 
baryon for non-uniform matter, nucleus (solid line) and shell 
state (dot-dashed line), and uniform matter (dashed line) are 
shown as functions of density for 1^ = 0.5 and T — 1 MeV. 
Lower panel (b): free energy difference per baryon between 
nucleus and shell states versus baryon density. 



uniform nuclear matter over a large range of baryon den- 
sities and 'with proton fractions, Yp = 0.2, 0.3, 0.4 and 
0.5. Within the spherical Wigner-Seitz approximation 
for non-uniform nuclear matter, we studied the system 
at a lo-w temperature, T = 1 MeV. We find ne'w shell 
states 'which minimize the free energy over a significant 
range of baryon densities. These shell states minimize the 
Coulomb energy of configurations ■with many protons at 
the expense of a larger surface energy. These states are 
related to a possible central depression in the density dis- 
tribution of super heavy nuclei. As the density increases, 
the system undergoes phase transitions from normal nu- 
clei to shell states, and then to uniform matter. 
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FIG. 10; (Color online) Free energy per baryon versus 
Wigner-Seitz cell radii for baryon densities ol ps — 0.040, 
0.043, 0.060, and 0.084 fm~^, proton fraction Yp = 0.3 and T 
= 1 MeV. The two open circles denote minima corresponding 
to nucleus (left) and shell (right) states. 
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FIG. 12; (Color online) The upper panel (a) shows atomic 
(solid line) and proton (dashed line) numbers, while the lower 
panel (b) shows WS cell radii, for nucleus (squares) and shell 
states (filled circles) as functions of the baryon density for Yp 
= 0.3. 
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FIG. 11: (Color online) The neutron (solid line) and proton 
(dash line) density distributions for the Wigner-Seitz cells at 
baryon densities of ps = 0.040, 0.043, 0.060, and 0.084 fm"^, 
for proton fraction Yp — 0.3 and T = 1 MeV. The red curves 
give density distributions for shell states. The black curves 
give density distributions of nucleus states. The correspond- 
ing cell radii are shown in each panel. The dotted lines in the 
lower right panel are the neutron (upper) and proton (lower) 
density of uniform matter. 



We compared our Wigner-Seitz results to two stabil- 
ity analyses of uniform nuclear matter. Here we looked 
at collective modes in the RMF plus Vlasov formalism 
and in the finite temperature RPA. The two analyses 
agree well with each other and give a slightly lower tran- 
sition density to uniform matter at large proton fractions 
{Yp = 0.5, 0.4), and a higher transition density at small 
proton fractions (below 0.3) compared to the WS results. 
This latter difference may suggest the existence of non- 
spherical configurations in the density region between 



FIG. 13; (Color online) Upper panel (a); the free energy per 
baryon for non-uniform matter, nucleus (solid line) and shell 
state (dot-dashed line) , and uniform matter (dashed line) are 
shown as functions of density for Yp — 0.3 and T — 1 MeV. 
Lower panel (b); free energy difference per baryon between 
nucleus and shell states versus baryon density. 

the spherical WS transition density and the collective 
modes/RPA transition density for non- uniform neutron 
rich matter. 

We conclude that even in a spherical approximation, 
there are more possible configurations than have been 
previously considered. These new shell states may play 
a role in the thermal and quantum fluctuations of the 
system and in its response to a variety of probes. Fur- 
thermore these shell states may allow one to include, in 
a simple way, much of the physics of, potentially compli- 
cated, non-spherical pasta configurations. 

In this investigation pairing between nucleons is not in- 
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FIG. 14: (Color online) Free energy per baryon versus 
Wigner-Seitz cell radius for baryon densities of pB ~ 0.035, 
0.044, 0.055, and 0.071 fm~^, proton fraction Yp = 0.2 and 
T = 1 MeV. The two open circles denote minima of nucleus 
(left) and shell (right) states. 



a more complete range of baryon densities, temperatures, 




FIG. 16: (Color online) The upper panel (a) shows atomic 
(solid line) and proton (dashed line) numbers, while the lower 
panel (b) shows WS cell radii, for nucleus (squares) and shell 
states (filled circles) as functions of the baryon density for Yp 
= 0.2. 
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FIG. 15: (Color online) The neutron (solid line) and proton 
(dashed line) density distributions for the Wigner-Seitz cells 
at baryon densities of ps = 0.035, 0.044, 0.055 and 0.071 
fm"^, with proton fraction Yp = 0.2 and T = 1 MeV. The red 
curves give density distributions of shell states. The black 
curves give density distributions of nuclear states. The corre- 
sponding cell radii are shown in each panel. The dotted lines 
in the lower right panel are the neutron (upper) and proton 
(lower) densities of uniform matter. 



eluded, since it is expected to have only a small effect on 
bulk properties such as the pressure or free energy. How- 
ever the effects of pairing on transport properties remains 
to be investigated. In the future, we will study the equa- 
tion of state at higher temperatures and investigate the 
response of shell states to supernova neutrinos. We will 
also study the system at smaller proton fractions appro- 
priate for neutron star matter in beta equilibrium. Our 
ultimate goal is to produce an equation of state spanning 
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FIG. 17: (Color online) Upper panel (a): the free energy per 
baryon for non-uniform matter, nucleus (solid line) and shell 
state (dot-dashed line) , and uniform matter (dashed line) are 
shown as function of density with Yp = 0.2 and T = 1 MeV. 
Lower panel (b): free energy difference per baryon between 
the nucleus and shell states versus baryon density. 



and proton fractions, that will be suitable for astrophys- 
ical applications to neutron stars and supernovae. 
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FIG. 18: (Color online) The phase diagram for nuclear mat- 
ter at T = 1 MeV and different proton fractions, combining 
results from Wigner-Seitz cell calculation, collective modes 
analysis and RPA at finite temperature. 
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